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ABSTRACT 

Several  boundary  integral  techniques  are  available  for  the  computation  of 
the  solution  to  Laplace's  equation  in  multi-connected  domains.  However,  for 
cases  where  the  domain  is  changing,  such  as  in  incompressible,  inviscid  fluid 
flow  with  free  surfaces,  iterative  methods  are  highly  attractive.  The  paper 
describes  one  such  formulation  and  tests  it  on  circular  and  elliptic  annuli. 

It  is  necessary  to  use  interpolated  quadrature  points  to  maintain  accuracy 

A 

when  regions  of  the  annuli  are  thin. 
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SIGNIFICANCE  AND  EXPLANATION 

Several  boundary  Integral  techniques  are  available  for  the  computation  of 
the  solution  to  Laplace's  equation  in  multi-connected  domains.  However,  for 
cases  where  the  domain  is  changing,  such  as  in  incompressible,  invlsctd  fluid 
flow  with  free  surfaces,  iterative  methods  are  highly  attractive.  Hie  paper 
describes  one  such  formulation  and  tests  it  on  circular  and  elliptic  annuli. 

It  is  necessary  to  use  interpolated  quadrature  points  to  maintain  accuracy 
when  regions  of  the  annuli  are  thin. 
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BOUNDARY  INTEGRAL  TECHNIQUES  FOR  MULTI-CONNECTED  DOMAINS 


G.  R.  Baker*  and  M.  J.  Shelley* 


Section  I.  Introduction 

Various  numerical  techniques  are  available  to  compute  solutions  to 
elliptic  partial  differential  equations.  For  specific  equations,  such  as 
Laplace's  equation,  the  biharmonic  equation  and  Helmholtz's  equation,  boundary 
integral  techniques  have  several  advantages  over  standard  finite-difference 
and  finite-element  techniques.  Highly  accurate  solutions  for  even  severely 
deformed  geometries  can  be  obtained  easily  by  boundary  integral  techniques. 

For  exterior  problems,  the  far-field  asymptotic  boundary  conditions  are 
automatically  satisfied.  No  special  effort  Is  required  when  the  domain 
changes  in  time  (or  with  some  other  independent  parameter),  since  points  on 
the  boundary  can  be  advanced  in  a  straightforward  fashion.  In  particular, 
generalized  vortex  methods,  based  on  boundary  integral  techniques,  have  been 
used  successfully  to  compute  free  surface  motion  of  inviscid,  incompressible 
fluid  (Baker  et  al.  1980,  1982). 

Several  other  researchers,  such  as  Longuet-Hi ggins  and  Cokelet  (1976)  and 
Pullin  (1982),  have  also  used  boundary  integral  techniques  to  study  free 
surface  motion.  They  used  costly  direct  matrix  inversion  techniques  to  solve 

3 

the  integral  equations,  which  take  0(N  )  operatons  to  perform  where  N  is 
the  number  of  points  that  represent  the  boundary.  In  contrast.  Baker  et  al. 
(1982)  realized  that  a  suitable  choice  of  source  or  dipole  distributions  along 
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the  surface  will  lead  to  Fredholm  integral  equations  of  the  second  kind  that 
may  be  solved  iteratively  in  0(N  )  operations.  In  simply-connected  domains, 
the  solution  to  Laplace's  equation  with  Dirichlet  boundary  conditions  may  be 
found  iteratively  when  a  dipole  distribution  along  the  boundary  is  used.  In 
multi-connected  domains,  an  external  source  contribution  must  be  added  to  the 
dipole  distribution  along  the  boundary  as  a  representation  for  the  velocity 
potential.  This  paper  will  describe  how  iterative  techniques  may  be  used  to 
find  both  the  source  strength  and  the  dipole  distribution.  Several  test 
examples  are  solved  numerically.  In  general,  standard  quadrature  techniques 
provide  accurate  solutions  to  the  boundary  integral  equations.  However,  when 
the  multi -connected  domain  involves  two  non-intersecting  surfaces  that  lie 
close  together,  interpolated  quadrature  techniques  are  used  to  solve  the 
boundary  integral  equations  accurately. 

One  application  of  these  results  is  to  the  study  of  accelerating  thin 
fluid  shells.  Baker  (1983)  has  shown  that  the  motion  may  be  determined  from 
the  solution  of  Laplace's  equation  in  a  multi -connected  domain  with  Dirichlet 
boundary  conditions.  Here  the  presence  of  the  source  term  is  crucial  in 
determining  correctly  the  dynamics  of  the  motion.  In  some  cases,  it  is  known 
on  other  grounds  that  there  is  no  flux  across  the  surfaces  and  so  no  source 
term  is  necessary.  This  is  relevant  in  studies  of  the  classical  Rayleigh- 
Taylor  instability  (Baker  et  al.,  1980)  and  of  the  motion  of  water  waves  over 
variable  bottom  topography  (Baker  et  al.,  1982). 

In  Section  II,  we  start  by  considering  Laplace's  equation  exterior  to  a 
simply-connected  domain.  This  simple  case  contains  the  essential  features  of 
the  application  of  boundary  integral  equations  for  potential  problems  in 
multi -connected  domains.  In  Section  III,  we  consider  Laplace's  equation 
between  two  closed  non-intersecting  surfaces  and  finally,  in  Section  IV,  a 
numerical  technique  and  results  are  presented. 
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Section  II.  The  Exterior  Problem 

Consider  the  region  F  exterior  to  a  simply-connected  domain  D  with 
boundary  3D.  For  convenience,  the  boundary  will  be  assumed  to  have  a 
continuous  normal.  Boundaries  with  corners  introduce  a  slight  modification  in 
the  mathematics  (for  details,  see  Jaswon  and  Symm,  1977).  Suppose  that  the 
potential,  <j>,  which  satisfies  Laplace's  equation  in  F,  is  specified 
along  3D.  In  addition,  the  behavior  of  41  far  from  3D  must  be  specified.  In 
three-dimensions,  a  unique  solution  is  determined  by  requiring  that  41  decays 
algebraically  far  from  3D.  However,  in  two-dimensions  the  situation  is  more 
complicated.  A  unique  solution  can  be  found  if  4>  is  logarithmic  or  tends  to 
a  constant  at  infinity  but  not  both  (Kellogg,  1929). 

For  simplicity  of  presentation,  we  shall  assume  that  the  potential  is 
logarithmic  at  infinity  if  the  geometry  is  two-dimensional.  The  other  case  is 
treated  similarly.  According  to  classical  potential  theory,  4.  may  then  be 
expressed  in  terms  of  a  source  distribution  0  along  3D; 

♦  (p)  *  /  o(q)g(p,q)dq  ,  P  c  DuaD  (2.1) 

3D 

where  p  and  q  are  field  points  and  g(p,q)  is  the  Green's  function  for 
Laplace's  equation  in  free-space.  In  particular,  when  p  e  3D,  equation  (2.1) 
constitutes  a  Fredholm  equation  of  the  first  kind  for  the  source  strength  a  in 
terms  of  the  specified  potential  4.  If  collocation  and  numerical  quadrature 
are  used  to  solve  equation  (2.1),  a  matrix  equation  results  which  is  usually 
solved  by  direct  inversion  techniques.  We  are  not  aware  of  any  matrix 
splitting  that  leads  to  an  iteration  scheme  that  converges  globally,  that  is, 
converges  for  any  choice  of  30.  Once  a  is  determined,  4  can  be  evaluated 
in  F  via  (2.1). 
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Alternatively,  as  suggested  by  the  form  of  a  solution  for  an  interior 
Dirichlet  problem,  we  may  seek  to  express  4  in  terms  of  a  dipole 
distribution  u  along  3D.  However,  such  a  representation  is  not 
sufficient.  A  dipole  distribution  decays  to  zero  at  infinity,  and  a  source 
term  must  be  added  in  order  to  satisfy  the  condition  there.  The  location  of 
the  source  inside  3D  is  arbitrary  but  for  numerical  computations  it  is  best 
not  to  place  it  near  3D.  For  convenience,  the  source  may  be  considered  to  be 
at  the  origin  of  the  coordinate  system.  Thus 

*(P)  *  /  u(q)  (p.q)dq  +  A<|>  (p)  ,  p  e  UusD  (2.2) 

where  n  is  the  normal  pointing  into  U  and  the  normal  derivative  is  taken 
with  respect  to  q  (hence  the  q  subscript  on  n).  The  unit  source 
potential  $s(p)  depends  on  the  nature  of  3D  and  the  spatial  dimension.  In 
particular,  for  a  closed  3D  in  two  dimensions, 

♦s(p)  --5!  l°9|p|.  (2.3) 

Clearly,  equations  for  v  and  A  must  be  sought. 

As  p  approaches  3D  along  the  normal,  (2.2)  takes  on  the  limiting  form 

I f  u(q)  42-  (p»q)dq  -  =  ${p)  -  A.t>  (p)  =  R  (p ) ,  p  c  3D  (2.4) 

3D  3  q  c  s 

Equation  (2.4)  has  been  written  explicitly  in  the  form  of  a  Fredholm  integral 
equation  of  the  second  kind  for  y.  The  arbitrariness  in  A  is  only  apparent 
in  that  the  Fredholm  alternative  must  be  satisfied  in  order  for  (2.4)  have  a 
solution  y.  As  required  by  Fredholm  theory,  the  eigenvalues  x  of  the 


equation. 


5 


2a  $  u (q )  (p,q)dq  -  p(p)  =  0»  (2.5) 

3D  °nq 

must  be  known.  Kellogg  (1929)  proves  that  the  eigenvalues  are  distinct  and 
real  and  they  all  satisfy  |a|  >  1.  In  addition,  a  =  1  is  an  eigenvalue  of 
geometric  multiplicity  one  with  eigenvector  y  =  C,  a  constant  (Jaswon  and 
Symm,  1977).  This  corresponds  to  a  potential  distribution  $  -  C  in  D  and 
<j>  =  0  in  TT. 

Since  (2.4)  has  a  non-trivial  solution  when  R  =  0  there  is  no  solution 
to  (2.4)  for  R  *  0  unless  P  satisfies  the  Fredholm  alternative.  Multiply 
(2.4)  by  a  source  distribution  a(p)  along  3D  and  integrate  around  30  with 
respect  to  p.  The  result  may  be  written  as 

/  u(q)[f  o(p>!jMp,q)dp  -  ^ijdq  *  /  R(pMp)dp  (Z.6) 

3D  3D  3nq  3D 

In  particular,  if  a  is  a  nontrivial  solution  of  the  adjoint  problem, 

r  o(q)  $Mp.q)dq  -  °.  (2,7) 

3D  ,np 

then 

/  R(p)o(p)dp  =  0,  (2.8) 

3D 

where  the  following  relationship  has  been  used; 

g(p.q)  *  g(q,p).  (2.9) 

Clearly,  (2.8)  is  a  necessary  condition  for  y  to  exist.  The  fact  that  it  is 
also  a  sufficient  condition  follows  from  the  usual  method  of  proof  that 
establishes  the  Fredholm  alternative  (see  Mikhlin,  1957).  The  Fredholm 
alternative  also  guarantees  that  there  exists  only  one  nontrivial  a  which 
satisfies  (2.7)  aside  from  a  multipl icati ve  constant. 
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Provided  R  satisfies  (2.8),  there  is  a  solution  y  for  (2.4)  which  is 
arbitrary  to  the  addition  of  any  constant  y.  Fortunately,  a  constant  y~ 
corresponds  to  a  constant  potential  in  D  and  zero  potential  in  IT.  In  most 
physical  applications,  such  as  incompressible,  inviscid  irrotational  fluid 
flow,  it  is  v<j>  that  is  important  and  any  constant  jT  may  be  ignored. 

Finally,  note  that  condition  (2.8)  becomes,  upon  substitution  for  R 
from  (2.4), 


A  /  4>c  (p )cr(p )dp  =  /  ♦(pMp)dp  (2.10) 

3D  s  3D 

which  is  an  equation  for  A.  Thus,  a  nontrivial  a  is  determined  from  (2.7) 
and  (2.10)  is  then  used  to  determine  A.  Next  (2.4)  is  solved  for  y  with  any 
additional  constraint  that  eliminates  the  arbitrariness  in  y.  For  example,  one 
may  specify  the  value  of  y  at  a  point  on  3D; 

y(p0)  •  0.  (2.11) 

At  first  sight,  equations  (2.4),  (2.7)  and  (2.10)  appear  to  introduce 
greater  computational  complexity  than  (2.1).  However  both  (2.4)  and  (2.7)  may 
be  solved  iteratively.  Let  y^(p)  be  the  nth  iterative  and  obtain  the  next 
iterative  from 


u(n+1)(p)  =  T.i(n)(p)  -  Ty(n)(p0)  (2.12a) 


where 

Ty(p)  =  2f3Qu(q)|^-(p,q)dq  -  2R(p). 


(2.12b) 
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As  n  +  «,  •*  y,  the  solution  to  (2.4)  subject  to  (2.11).  Similarly,  an 

iteration  procedure  for  a  may  be  used  where  for  convenience  the  additional 
constraint, 

max  I  a (p ) I  =  1  (2.13) 

pe3D 

is  imposed  to  remove  the  arbitrary  multiplicative  constant  in  a.  The 
convergence  of  the  iteration  scheme  follows  from  the  global  convergence  of  the 
Neumann  series  for  both  integral  equations.  A  proof  follows  from  the  proof 
given  in  Baker  et  al  (1982)  for  the  case  of  open  periodic  surfaces.  In 
particular,  for  time  dependent  domains  the  iteration  is  very  efficient  since  a 
good  first  guess  is  always  available  from  information  at  the  previous  time 
value,  and,  if  information  from  previous  time  levels  is  retained, 
extrapolation  further  improves  the  first  guess.  If  the  solution  <j>  to  the 
two-dimensional  Laplace's  equation  is  being  sought  in  D  subject  to  <j>  tending 
to  a  constant  at  infinity,  4>s (p )  =  1  must  be  used  in  place  of  (2.3)  and  the 
source  distribution  in  (2.1)  must  be  modified  (Jaswon  and  Symm,  1977). 

As  shown  in  the  simple  case  above,  the  Fredholm  alternative  lies  at  the 
heart  of  the  application  of  dipole  distributions  to  the  solution  of  Laplace's 
equation  in  multi -connected  domains.  In  the  next  section,  a  more  general 
multi -connected  domain  will  be  considered  and  then  numerical  results  will  be 
presented  in  the  following  section. 

Section  III  The  Annular  Problem 

Consider  the  multi -connected  domain  as  shown  schematically  in  Figure  1. 
Domain  02  lies  between  two  closed,  non-intersecting  surfaces 
30^  and  aD,,  with  aD^  enclosing  aO^.  The  exterior  domain  D  is  composed  of 
two  parts,  U3  internal  to  aDg  and  Dj  external  to  aD^.  Once  again,  aDj 
and  3D2  are  assumed  to  have  continuous  normals. 


■  I  "V  «  v  ' 


■  .'v  '.■*  .V  .VA-.T  ■:■  V  ..■  ^ 


.  I". 
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To  be  specific  and  for  convenience,  suppose  that  a  solution  to  the  two 
dimensional  Laplace's  equation  is  sought  in  D2  with  Dirichlet  boundary 
conditions  imposed  at  30^  and  303.  Following  the  procedure  adopted  in 
Section  II,  the  solution  is  expressed  in  terms  of  a  source  tern  and  dipole 
distributions  along  30^  and  3D2; 

2 

<t*(p)  ■  -^logjpj  +  ^  /  Uj  (q )  |^~(p,q)dq,  p  e  0?  .  (3.1) 

J 

Let  3Dj  and  302  be  parameteri sed  in  a  counter-clockwise  direction  by 

(xi(e).yi(e))  and  (x2(e),y2(e))  respectively.  In  free-space,  the  two- 

1  2  2 

dimensional  Green's  function  i s  1 og{ (x-x^  )  +  (y-yk)  >  where  (x^.y^) 

locates  the  source  point  on  3Dk  .  The  normal  derivative  of  the  Green's 
function  evaluated  at  the  kth  surface  has  the  form 


„  1  xke  '  Hyj  (e  )-yk  (e  '  )) _yke  (e  '  (e)-xk (e ' )} 

“  -n~~  - t - - n - 

J  {Xj  (e)-xk  (e '  )r+{yj  (e)-yk(e  ' )} 


(3.2) 


where  the  subscript  e  denotes  differentiation  and  the  field  point  lies  on 
the  jth  surface. 

The  evaluation  of  (3.1)  at  3D^  and  3D2  gives  two  coupled  Fredholm 
integral  equations  for  u ^  and  V2: 


X  f  ui (e' )K,. (e,e' )de '  +  X  /  u?(e ' )K  (e, e ' )de '  + 
3Dj  3  D  g 

=  ^i(e)  -  log  (x^(e)  +  y^(e)} 


(3.3a) 


interpolated  points 
^original  points 
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defined  by  the  requirement  that  (x^eg),  yk(eg))  is  the  closest  point  on 
surface  k  to  (xj(e),  yj(e)).  Figure  3  gives  a  schematic  of  the 
situation.  The  value  eg  may  be  found  by  using  interpolation  and  Newton 
iteration  to  locate  the  minimum  distance.  To  resolve  the  peak  or  trough,  N 
points  are  interpolated  along  the  surface  in  the  following  v/=>y.  Given  that 
points  (xk(e‘),  yk(e'))  are  originally  assigned  by  evenly  spaced  intervals 

A 

in  e’,  N  new  points  are  obtained  by  taking  evenly  spaced  intervals  in  e, 
where 

e'  =  e  +  a  sin(e  -  e'g)  (4.13) 

A 

and  e  =  e'g  gives  one  of  the  new  points.  The  trapezoidal  rule  is  used  with 
the  new  points;  quantities  such  as  t  and  y  are  also  interpolated  to  be 
able  to  evaluate  the  integrand.  Note  that  the  same  number  of  points  N  are 
used.  Clearly,  the  choice  of  a  will  dictate  the  accuracy.  Note  that 
interpolation  must  be  done  whenever  a  point  on  the  jth  surface  is  too  close 
to  the  kth  surface. 

The  error  E  for  the  circular  case  may  be  analysed  in  more  general  terms 

to  estimate  its  behavior  for  other  geometries.  For  large  N,  the  error  T  i 

-N 

most  strongly  affected  by  the  term  p  ,  provided  p  t  1.  When  j  -  1  and 
k  =  2,  the  error  is  associated  with  the  approximation  to  the  integral 
involving  K^e.e').  Set  =  R  and  R|  =  R+D.  Then, 

p'N  =  (1  +  D/R)"N 

=  exp(-N  1 og ( 1  +  D/R))  (4.14) 


exp  (-2ttD/aS  (1  +  D/2R  ) ) 
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is  quite  small.  Incidentally,  the  mode  m  =  16  cannot  be  resolved  by  the 
numerical  quadrature  since  alternate  points  were  used  so  m  =  15  is  the 
highest  mode  realistically  treated  by  the  quadrature.  The  maximum  number  of 
iterations  required  occurs  when  m  =  1  and  the  number  decreases  for  larger 
values  of  m,  consistent  with  other  iterative,  numerical  techniques  used  to 
solve  elliptic  problems.  Clearly,  multi-grid  techniques  may  be  used  to  reduce 
strongly  the  required  number  of  iterations,  but  we  have  not  pursued  that 
aspect  here.  The  numerically  calculated  value  for  the  source  strength  A  was 
always  within  roundoff  error  of  its  exact  value. 

Unfortunately  there  are  times  when  straightforward  use  of  the  numerical 
procedure  described  above  will  result  in  errors  of  0(1).  From  (4.12),  it  is 
easy  to  see  that  E  is  0(1)  when  Rj  and  R2  are  very  close  together. 
Figure  2  displays  the  behavior  of  E  as  the  width  of  the  circular  annulus  is 
decreased  in  the  test  case  (4.4-6).  Here  Rj  =  1+D,  R^  =  1,  m  «  1,  Nj  =  N2  = 
N,  and  an  absolute  tolerance  of  10"^  was  used  for  the  convergence  of  the 
iterations,  (3.7)  and  (3.8).  Clearly  for  a  fixed  number  of  discretisation 
points  N,  the  strong  variation  in  the  integrand  Kjk,  j  ^  k,  is  not  well 
resolved  when  a  point  on  one  surface  is  close  to  the  other  surface  along  which 
the  integral  is  performed.  The  error  begins  to  downgrade  noticeably  when  the 
thickness  D  of  the  annulus  is  sufficiently  small.  This  error  has  been 
previously  observed  by  Maskew  (1977)  in  a  related  calculation  involving  vortex 
sheet  motion.  Clearly,  as  Maskew  points  out,  more  points  are  required  to 
resolve  the  strong  variation  of  the  integrand,  and  as  N  is  increased  for 
fixed  D  accuracy  improves  dramatically. 

It  is  expensive  and  unnecessary  to  increase  N  substantially  for 
accuracy  when  D  is  small  .  Consider  a  point  (xj (e),yj (e))  near  surface 
k.  The  integrand  Kjk(e,e')  has  a  sharp  peak  or  trough  centered  around  e^  , 


N 

m 

Error 

I 

16 

1 

0.305  x  10“3 

26 

32 

1 

0.466  x  10-8 

26 

64 

1 

rH 

f-H 

1 

o 

rH 

X 

r*» 

H 

rH 

• 

o 

26 

32 

2 

0.528  x  10-8 

14 

32 

4 

0.160  x  10-7 

8 

32 

8 

0.239  x  10-6 

5 

32 

12 

0.382  x  10'5 

4 

32 

15 

0.305  x  10"4 

3 

Table  I:  Results  for  the  circular  annulus  as  described  in  the  test 
in  all  cases. 
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w2  *  (2  +  F)  (cos  (jmh)  -  1)  (4.11b) 


where 


E 


o/  m  ,  ~rn « 

2(p  +  P  ) 


N  N-fli  m 
P  -  P  -  P  - 


(4.12) 


and  p  *  R ^ /R 2  . 

A  numerical  code  has  been  written  that  solves  the  Dlrlchlet  problem  for  a 
general  annulus  (not  necessarily  circular).  In  particular,  the  code  was 
applied  to  the  test  case  described  above.  Results  are  shown  in  Table  I  when 
R1  *  2»  R2  =  and  N1  *  N2  =  N  points  are  used.  The  reported  error  E  is 
the  maximum  absolute  difference  between  the  computed  dipole  strengths  and  the 
exact  values  given  In  (4.7);  theoretically  this  error  should  be  2|F|  and  the 
agreement  Is  perfect  as  long  as  the  error  is  above  the  absolute  tolerance  of 
10“12  used  to  determine  convergence  of  the  iteration  scheme.  Also  tabulated 
are  the  number  of  Iterations  1^  and  Iy  required  to  solve  the  integral 
equations  for  and  pj  respectively  to  within  the  absolute  tolerance.  We 
emphasize  that  the  error  arises  solely  from  the  numerical  approximation  to  the 
integral  that  determines  the  contribution  to  the  field  point  on  one  surface 
from  the  dipole  distribution  along  the  other  surface.  For  m  =  1,  there  is  a 
dramatic  improvement  in  accuracy  as  N  increases  reflecting  the  infinite 
order  of  the  trapezoidal  rule  on  periodic  integrands.  Already,  for  N  =  64, 
the  error  is  below  the  tolerance  required.  Of  course,  this  result  is  not  too 
surprising  since  a  circular  annulus  has  been  chosen,  but  later  results  will  be 
represented  for  elliptic  annuli  that  still  show  the  high  accuracy  of  the 
method.  As  the  mode  number  m  of  the  potential  is  increased  with  fixed  N, 
the  resolution  of  the  integrand  deteriorates,  but  even  for  m  =  15,  the  error 
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The  eigenvectors  for  the  discrete  equations  that  represent  (3.5)  are 
proportional  to  cos(mjh)  and  sin(mjh).  In  particular,  for  m  =  0,  the 
eigenvalue,  which  corresponds  to  x  =  1  in  the  continuous  case,  is 


X  = 


R;  -  rn2 


(4.9) 


and  the  numerically  calculated  eigenvector  is  tj  =  1,  t2  =  -1,  which  is 
exact.  When  these  values  for  t^,  t2  are  substituted  into  (3.4)  and  the 
trapezoidal  rule  is  used  to  compute  the  integrals  in  (3.4),  the  numerically 
determined  value  for  A  is  also  exact,  that  is  A  =  1. 

Next,  the  error  in  solving  (3.3)  can  be  computed.  The  following  sums. 


Re 


N-l 


k+j  odd 


cos(kmh)  eikh  , 

i  jh  ’ ikh  i 
e  -  e 


for  m  *  0 
for  m  =  0 


(4.10a) 


Re 


Ny*  cos  (kmh)  eil(h  , 
k=0  Rje’J"  -  R^e 11(11 


1 

1 


m  r,N-m 


R'  R 


1  2 


RN-m  Rm 
K1  K2 


cos(jmh),  m  t  0 

(4.10b) 


give  the  trapezoidal  approximation  for  the  various  integrals  in  (3.3).  Once 
again,  (4.10a)  implies  that  the  principal -value  integrals  are  computed  exactly 
and  errors  come  only  from  the  approximation  to  the  other  integrals.  It  is  a 
straightforward  calculation  to  find  the  numerically  determined  dipole  sheet 
strengths, 


uj  =  (2  +  E)  (cos  (jmh)  -  1) 


(4.11a) 
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^(e)  =  (1  -  (|q-)  )  cos(me)  +  log  (Rj)  (4.6a) 

R  ^ 

*2(e)  =  ((j^)  -  1)  cos  (me)  +  log  (R2)  (4.6b) 

on  the  boundaries  30^  and  3D2  respectively.  The  exact  values  for  the 
dipole  sheet  strengths  such  that  y^(o)  s  p2(o)  =  0  are 

uj(e)  *  y2(e)  =  2  Ccos(me)  -  1].  (4.7) 

For  this  simple  test  case,  the  error  Involved  in  using  trapezoidal 
quadrature  can  be  calculated  exactly.  The  following  sums 


9  i  ,-h  N-i 
Re  {  £  e1Jh  l 
k=0 


Re  {  Rieijh  l 


k+j  odd 
N-l 


cos  (kmh) 


cos(kmh) 


}  •  { 


for  m  ^  0 
for  m  *  0 


(4.8a) 


R%N-m+  R^" 


k=0  R1e1Jh-R2e11 


1  12 


D  '»  n  ” 
K1  '  K2 


cos(jnh),  m  #  0 


(4.8b) 


where  h  =  2n/N  and  0  <  m  <  ,  give  the  trapezoidal  approximation  to  the 

various  Integrals  in  (3.5)  where  N  evenly  spaced  points  have  been  used  to 
represent  each  surface.  The  sum  (4.8a)  which  approximates  the  principal -value 
Integrals,  gives  the  exact  analytic  result,  reflecting  the  spectral  accuracy 
of  the  trapezoidal  rule.  Thus  the  only  errors  arise  from  the  sums  (4.8b)  that 
approximate  the  integrals  along  the  surface  opposite  to  the  one  containing  the 
field  point,  (xj(e),  y^?)),  in  K jk (e, e ' )  in  (3.2). 
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The  trapezoidal  rule  is  applied  as  follows.  Nj  points  are  selected  at  evenly 
spaced  intervals  in  e  to  represent  surface  3D..  At  each  such  point,  e^ 

*j 

say,  in  (3.7)  and  (3.8)  as  modified  by  (4.2)  and  (4.3),  Tj (e^ )  must  be 
calculated.  One  integral  involves  contributions  from  the  other  surface  and 


the  trapezoidal  rule  may  be  applied  directly.  For  the  integral  along  3D- 
however,  application  of  the  trapezoidal  rule  would  require  evaluation  of  the 
integrand  at  e‘  =  ej  which  is  an  indeterminant  form.  The  limiting  value  is 
easily  calculated  but  involves  derivatives  of  the  dipole  sheet  strength  and 


the  surface.  Instead,  it  is  more  convenient  to  apply  the  trapezoidal  rule  on 
alternate  points  so  that  e‘  =  ek  for  all  k  such  that  k  +  j  is  odd. 


Clearly,  a  quadrature  point  never  falls  on  e  and  the  integrand  is  always 
easily  computed.  A  disadvantage  to  this  procedure  is  a  slight  loss  of 
resolution  of  the  integrand.  Tests  will  show  that  this  is  not  important  in 
most  cases. 

As  a  simple  test,  consider  the  annulus  lying  between  the  circles  defined 
by 


and 


Xj  =  Rj  cos  (e)  ,  y^  RjSin  (e)  (4.4a) 

x2  =  R2  cos  (e)  ,  y2  =  R2  sin  (e)  (4.4b) 


where  Rj  and  R^  are  constants.  A  solution  to  Laplace's  equation  in  polar 
coordinates  (r,  9  )  in  the  annulus  is 


R  171 

(-£)  )  cos  (me)  +  log(r) 


(4.5) 


which  takes  on  the  values 


I 


Thus  the  integrals  in  (3.7)  may  be  replaced  by 


and 


-2/o  {«/"’(<!  ')  -  v\n)  (e»  Kn(e,e')  de'  -  Ul(n,(e) 

-2/  u2(n)<e')  K12(e,e')  de'  (4.2a) 

3D  2 


2  /  {u2(n)(e')  -  w2(n)(e)>  Ue.e1)  de' 

3D2  c  cc 

*2  {  u,!n)(e')  K2,(e,e')  de' 

3DX 


+  uz(n)(e) 

(4.2b) 


These  integrals  have  smooth,  periodic  integrands  and  may  thus  be  evaluated 
accurately  by  the  trapezoidal  rule.  Indeed,  the  accuracy  is  infinite  order  or 
spectral  (Isaacson  and  Keller,  1966).  Similarly,  the  integrals  in  (3.8)  may 
be  replaced  by 


/  (Tj(n)(e')  Gn(e,e')  -  M(n)(e)  Ku(e.e')>  de’  +|x1(n,(e) 

3D  j 


+  /  T?(n)(e' )  G12(e,e' )  de' , 

3D2 

j  {T2(n)(e‘ )  622(6,6')  -  r2(n)(e)  K22(e,e')}  de*  -  }  x2(n)  (e) 

3^2 

-  /  Ti(n)(e')  G21(e,e')  de'  . 


(4.3a) 


(4.3b) 
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T2(Tl’T2^e)  =  “  /  Ti(e')  G21(e,e')  de'  - 
3D  ^ 

( 3.8d) 

?  (e ' )  Gpp(e,e‘ )  de'  -  T2(e) 

3D2  c 

and  T  =  max  T .(T^n^,Tpn^)(e).  The  disadvantage  to  shifting  the  eigenvalues 
m  j=1>2  J  i  ^ 

is  that  the  convergence  rate  is  halved. 

The  above  interation  method  is  effectively  a  Neumann  iteration  method 
applied  to  the  discrete  version  of  the  boundary  integral  equations.  A  referee 
has  pointed  out  that  a  generalized  conjugate  residual  algorithm  might 
significantly  improve  the  convergence  (see  Eisenstat  et  al.  1983).  The  basic 
conjugate  gradient  method  was  tried  on  the  exterior  problem  with  little  gain 
in  convergence,  but  we  have  not  as  yet  tried  the  improved  algorithm. 

Clearly,  the  numerical  implementation  of  the  above  procedure  demands 
reliable  and  accurate  evaluation  of  the  integrals.  In  the  next  section,  a 
numerical  quadrature  and  results  for  some  test  cases  are  presented. 


Section  IV.  Numerical  Results 

The  design  of  a  numerical  quadrature  for  the  integrals  in  (3.7)  and  (3.8) 
depends  on  several  factors.  Firstly,  the  kernels  Kjj  and  Gj j  are 
singular.  Secondly,  the  smoothness  of  the  surfaces  and  and  of  the 

boundary  values  <|>j  and  <j>2  may  affect  the  accuracy.  Fortunately,  for 
applications  to  studies  of  free  surface  flows,  the  surfaces  and  boundary 
values  are  generally  C°°  functions.  Furthermore,  for  closed  surfaces  in  two- 
dimensions,  the  singularities  in  Kjj  and  Gjj  may  be  removed  by  using  the 
identity. 


L2 


‘  2*2(e)  +  Tt  109  {X2(e)  +  y2(e)} 


will  converge  to  values  of  and  ^2  that  are  shifted  by  some  constant 
value  from  the  correct  values.  However  the  net  effect  of  these  shifts  is  to- 
produce  a  corresponding  constant  shift  in  the  potential  distribution. 

Normally,  such  shifts  in  potential  are  not  physically  interesting.  However, 
if  necessary,  it  Is  a  relatively  simple  matter  to  add  the  appropriate  constant 
to  the  potential. 

The  iteration  procedure  for  finding  the  eigenfunctions  tj  and  t2  that 
correspond  to  x  -  1  must  be  modified  differently  since  the  eigenfunctions 
tj  and  t 2  that  correspond  to  X  =  -1  are  not  known  in  general.  Instead  the 
eigenvalue  x  =  -1  is  shifted  to  the  origin  by  using  the  following  iteration 
scheme ; 


T^(n  +  l)(e)  =  Tl(T{n),T^nh(e)/Tm  (3.8a) 

t  2  (n  +  ^(e)  =  T2(Tin^,T^n^)(e)/Tm  (3.8b) 


where 


Ti(Ti»T2)(e)  =  f  Ti^e')  G^(e,e')  de1  + 

8  L)  ^ 

(3.8c) 

/  T2(e')  G^2^e,e de'  +  Tj(e) 
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Gjk (e,e' )  =  -  K|<j(e'  e)  (3.6) 

and  t  is  related  to  a  source  distribution  o  by  t  =  osg  where  s  is  the 
arclength. 

The  solution  4  may  now  be  calculated  as  follows.  The  solution  to  (3.5) 
yields  and  t2  which  are  substituted  into  (3.4)  so  that  A  may  be 
calculated.  Thus  (3.3)  may  be  solved  for  y^  and  y2  and  4>  maY  be 
evaluated  from  (3.1). 

There  remains  the  question  of  how  to  compute  the  solutions  to  (3.3)  and 
(3.5).  The  observation  has  already  been  made  that  A  =  1  is  an  eigenvalue  to 
the  homogeneous  equations  associated  with  (3.3).  There  is  also  an 
eigenvalue  A  =  -1  with  ^  =  C,  y2  3  -C,  which  corresponds  to  a  potential 

distribution  <|>  =  0  in  Dj  and  D3  and  <j>  =  C  in  03.  Unfortunately  the 

eigenvalue  a  =  -1  affects  the  iteration  scheme  that  arises  naturally  by  the 
generalization  of  (2.12).  Since  all  other  eigenvalues  satisfy  |a|  >  1  (see 
Kellogg,  1929),  the  eigenvalues  A  =  +  1  are  the  only  ones  that  prevent 
convergence  of  the  iteration  scheme.  Thus  the  modified  iteration  scheme, 

u|n+1)(e)  =  T1(u{n),y^n))(e)  -  T^y}" } ,y£n } ) (0)  (3.7a) 

y£n+1)(e)  -  T2(y{n),y^n))(e)  -  T2(yjn),y <n))(0)  (3.7b) 


where 


T^(y^,U2)(e) 


-2f  y  1  (e '  )K^(e,e'  )de'  -2/  y2(e  '  )K12 (e,e  '  )de ' 
3Dj  902 


+  2^ x (e)  -  log  {x^(e)  +  y^(e)}. 


(3.7c) 


(3.3b) 


X  j  pj(e' )K21(e,e' )de'  +  X  £  P2(e' )K22(e,e' )de' 
3D  j  SDg 

=  "  37  (x2^e^  +  y2(e)) 


L. 


where  ^(e)  and  $2(e)  are  the  values  Imposed  on  $  at  30^  and  dD2 
respectively.  The  origin  of  the  coordinate  system  has  been  chosen  to  coincide 
with  the  location  of  the  source  point  inside  D3  of  strength  A.  Although 
X  =  1,  it  has  been  introduced  to  facilitate  the  discussion  on  the  existence  of 
solutions  for  p^  and  p2  and  the  global  convergence  of  the  Neumann  series  for 
the  integral  equations. 

When  x  s  1  ,  pi  *  0  and  p2  =  C,  a  constant,  are  solutions  to  the 
homogeneous  equations  ($,  ■  <j>2  =  A  =  0)  corresponding  to  a  potential 
distribution  4,  -  0  in  and  D2,  4  *>  C  in  Dg.  This  eigenvalue  is  closely 
associated  to  that  in  (2.5).  According  to  Fredholm's  theorem  of  the 
alternative,  (3.3)  has  a  solution  provided 

•jfe  [  /  T^(e)  log  {Xj(e)  +  y*(e)}  de  +  /  x2(e)  log  (x2(e)  +  y2(e)}de] 

3U^  3  D2 

=  /  ♦1(e)tj(e)  de  +  /  <1>2 (e )x2 (e )  de  ,  (3.4) 

3D^  3 

where  and  12  satisfy  the  homogeneous,  adjoint  equations 

T,(e) 

£  x^e')  Gjjfe.e')  de'  +  /  x2(e')  G12(e,e‘)  de' - 2 — =  0  (3.5a) 

dD  ^  3^2 

/  T,(e')  G21(e,e‘)  de'  +  f  x2(e')  G22(e,e')  de'  +  ■  -2-g— 


=  0  (3.5b) 


where  the  quantity  aS  is  the  arc  length  between  points  and  so  NaS  =  2nR. 

This  result  in  (4.14)  is  obtained  by  rewriting  D/R  =  2ttD/NaS  and  considering 
D/aS  fixed  for  large  N.  The  approximation  z/(l  +  z/2)  to  log  (1  +  z)  is 
used  rather  than  z ( 1  -  z/2).  We  found  that  this  modification  gave  better 
agreement  with  the  numerical  results  reported  later. 

For  the  other  case,  j  =2  and  k  =  1;  set  Rj  =  R,  R£  =  R-D  and 
NaS  =  2uR  Now, 

p“N  =  (1  -  D/R)+N  -  exp(-2irD/(AS(l  -  D/2R ) ) ).  (4.15) 

The  factors  1±D/2R  in  (4.14)  and  (4.15)  are  curvature  j<|  =  1/R 
corrections  to  the  expression  exp(-2irD/AS )  which  describes  the  influence  of 
the  error  as  a  field  point  approaches  a  flat  surface  containing  a  dipole 
distribution  (the  error  for  this  case  is  also  easily  calculated).  In  general, 
the  correction  can  be  written  as  1  +  D|<|/2,  1  -  D|k|/2  depending  upon 
whether  the  region  of  the  surface  nearest  the  point  appears  convex  (see  Figure 
2)  or  concave  respectively.  In  Figure  4,  we  replot  the  error  as  a  function  of 
D/(aS(1  -  D|k|/2))  where  as  and  k  is  determined  on  the  outer  surface. 
Clearly,  the  asymptotic  form  (4.15),  while  not  strictly  valid  for  all  values 
of  D,  does  capture  the  essential  behavior  of  E. 

To  decrease  the  error,  the  local  AS  between  the  new,  interpolated  points 


should  be  given  by 


aS(1-DM/2) 

Figure  4.  The  same  as  Figure  2  except  D  is  scaled  by  AS(l-d|ic|/2)  where 
the  local  arclength  between  points,  AS,  and  <  are  measured  on 
the  outer  boundary. 


where  C  is  a  constant,  chosen  by  experimentation  to  give  the  required 

accuracy.  From  (4.13),  it  is  easy  to  determine  the  new  AS  in  terms  of  a, 

AS  =  AS(1  +  o).  (4.17) 

Consequently,  a  has  the  form 

a  «  £^/(l  ±^)  -  1.  (4.18) 

In  Figure  5,  the  error  is  shown  as  D  is  varied  for  the  same  test  case 

as  shown  in  Figure  2,  but  now  interpolated  points  are  used  whenever  a  >  0  in 

(4.18).  Interpolation  was  done  using  a  discrete  Fourier  representation  for 
which  C  =  9.0  was  found  to  be  a  good  choice.  For  small  D,  the  errors  are 
essentially  the  same  as  for  larger  D  except  for  extremely  small  values  of 
D,  when  the  interpolated  points  are  all  clustered  near  e0  and  there  are 
effectively  no  points  resolving  the  rest  of  the  surface.  This  effect  clearly 
depends  on  N. 

In  some  cases,  such  as  a  "noisy"  representation  of  the  surface,  Fourier 
interpolation  is  inappropriate.  In  such  cases,  cubic  spline  interpolation  may 
be  used  and  we  show  the  results  in  Figure  6.  Because  of  the  errors  involved 

4 

in  using  cubic  splines,  the  error  can  be  no  better  than  0(Ae  ).  Now  there  is 
a  transition  from  spectral  accuracy  to  0(Ae4)  as  D  becomes  small  and 
interpolated  points  are  used.  For  this  case,  C  =  6.0  since  the  accuracy  is 
limited  and  there  is  no  point  in  using  interpolated  points  for  a  large  value 
of  C. 

Finally,  we  demonstrate  that  these  results  are  not  special  to  a  circular 
annulus  but  hold  true  for  generally  shaped  annuli.  Take  the  boundaries  of  the 


D  when  F< 


annuli  to  be  the  following  ellipses. 


Xj  =  r  cos(e) 

l/2 

yx  =  (r2  -  3/4)  sin(e). 


(4.19a) 

(4.19b) 


and 

(4.20a) 

(4.20b) 


X2  =  cos(e) 
y2  =  7  sin(e)  . 


Our  tests  keep  the  inner  boundary  fixed  and  vary  the  outer  boundary  by 
changing  r.  The  conformal  map. 


x  +  iy 


.  /3 
" 


cosh(c  +  in). 


(4.21) 


transforms  the  elliptic  annulus  into  a  circular  one  and  so  an  exact  solution 
for  the  Dirichlet  problem  can  be  found.  In  particular,  we  chose  the  boundary 
conditions. 


^  =  (a  cosh  (2^)  +  3  sinh(2c1))  cos  (2e)  (4.22a) 

4>2  =  (7  a  +  7  6)  cos  (2e)  ,  (4.22b) 


where 


30 


a  =  -6  -  u1/(cosh(251)  +  sinh^^)) 

(4.23a) 

0  4  - 

8  "  "lv2 

(4.23b) 

cosh(c. )  =  — 

1  /J 

(4.23c) 

sinhCcj)  =  (4r2/3  -  1)1/2 

(4.23d) 

and  the  corresponding  dipole  strengths  are 


nj  =  yj  cos(2e) 

(4.24a) 

U2  *  W2  cos(2e) 

(4.24b) 

First,  results  are  shown  in  Figure  7  for  the  elliptic  annulus  without  the 
use  of  interpolated  points  and  j7j  =  1  and  jTg  =  Here,  D  was  chosen  as 
the  minimum  thickness  of  the  annulus.  The  pattern  of  results  is  similar  to 
Figure  1.  Next,  we  replot  the  error  in  Figure  8  as  a  function  of 
D/(aS(1  -  D|k|/2)),  where  the  local  as  and  k  are  measured  on  the  outer 
surface,  and  find  that  the  error  behavior  is  again  well  described  by  (4.15). 
Finally,  in  Figures  9  and  10,  results  are  given  when  interpolated  points  are 
used,  Fourier  and  cubic  spline  interpolation  being  used  respectively. 

The  relative  merits  of  using  Fourier  and  cubic  spline  interpolation  hinge 
upon  accuracy  versus  cost.  Our  approach  is  to  use  a  fixed  number  of 
quadrature  points,  hence  the  added  cost  is  determined  by  the  type  of 
interpolation  used.  For  illustrative  purposes,  we  chose  two  disparate  types 
of  interpolation:  discrete  Fourier  and  periodic  cubic  splines.  While  the 
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aS(1-DM/2) 

Figure  8.  The  same  as  Figure  7  except  D  is  scaled  by  &S(1-d|k|/2)  where 
the  local  arclength  between  points,  AS,  and  k  are  measured  on 
the  outer  boundary  at  the  point  of  minimum  thickness. 
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Figure  9.  The  maximum  of  the  absolute  error  E  in  dipole  sheet  strength  for 

the  elliptic  annulus  as  a  function  of  D  when  Fourier  interpolation 
is  used  to  redistribute  quadrature  points. 


Fourier  interpolation  is  spectrally  accurate,  each  interpolation  involves  a 
summation,  and  the  cost  is  increased  by  a  factor  0(N).  Fortunately,  Fourier 
interpolation  is  fully  vectorizable,  considerably  reducing  the  cost.  When  the 
surface  representation  is  "noisy",  or  the  Fourier  interpolation  considered  too 
expensive,  periodic  cubic  splines  offer  an  attractive  alternative.  The  error 
is  now  at  best  0(Ae  ),  but  the  cost  is  only  increased  by  a  factor  0(1). 

Section  V.  Conclusion 

In  conclusion,  it  is  possible  to(solve  elliptic  problems  in  multi- 
connected  domains  with  smooth  boundaries  using  iterative  boundary  integral 
techniques.  High  accuracy  is  possible  for  relatively  few  boundary  points, 
even  when  the  multi  connected  region  has  thin  parts.  A  simple  redistribution 
of  points  insures  the  high  accuracy.  This  technique  is  now  being  applied  to 
the  study  of  accelerating,  thin  fluid  shells  (Baker  1983)  and  of  the  fluid 
motion  of  vortex  layers;  these  results  will  be  published  elsewhere. 
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